Stretched exponential behavior and random walks on diluted hypercubic lattices 
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Diffusion on a diluted hypercube has been proposed as a model for glassy relax- 
ation and is an example of the more general class of stochastic processes on graphs. 
In this article we determine numerically through large scale simulations the eigen- 
value spectra for this stochastic process and calculate explicitly the time evolution 
for the autocorrelation function and for the return probability, all at criticality, with 
hypercube dimensions iV up to iV = 28. We show that at long times both relax- 
ation functions can be described by stretched exponentials with exponent 1/3 and 
a characteristic relaxation time which grows exponentially with dimension N. The 
numerical eigenvalue spectra are consistent with analytic predictions for a generic 
sparse network model. 
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INTRODUCTION 



In 1854 R. Kohlrausch used a phenomenological expression 

q K (t) = eM-(t/rf) (1) 

to parametrize the non-exponential decay of the electric polarization of Leyden jars (primi- 
tive capacitors )D} his son F. Kohlrausch later used the same expression to analyse creep in 
galvanometer suspensions^. A century later, in 1951 Weibull introducecP the closely related 
Weibull function; this survival probability functiorP which is widely used in the engineering 
literature is strictly of the Kohlrausch form, Eqn. ([I]). In 1970 Williams and Watts re- 
discovered the Kohlrausch function in the context of dielectric relaxation^. Under the name 
of "stretched exponential"^ the KWW (Kohlrausch- Williams- Watts) function has become 
ubiquitous in phenomenological analyses of non-exponential relaxation data, experimental 
or numerical. In particular the KWW form was used by Ogielski in a phenomenological fit 
the decay of the autocorrelation function at equilibrium for a 3d Ising spin glass modeP. 

Many arguments have been given as to why under certain assumptions, specific systems 
should show KWW relaxatiorP^, but there have always been lingering suspicions that for 
most cases the KWW expression is nothing more than a convenient fitting function of no 
fundamental significance. 

It was conjectured^ that KWW relaxation is the signature of a complex configuration 
space. Thus from the argument which follows it was suggested that random walks on a 
diluted hypercube (a hypercube with a fraction p of vertices occupied at random) near the 
critical concentration for percolation would lead to an autocorrelation function decay 
of the form q(t) ~ exp[— (t/r)"], with a specific value of the exponent, (3 = 1/3. 

For random walks at percolation threshold in a randomly occupied Euclidean (flat) space 
of dimension d such as Z d , the familiar Fickian diffusion law (R 2 ) ~ t is replaced by a 
sub-linear diffusion (R 2 ) ~ t^ d , with /3d = 1/3 for d > (P^. Random walks on the surface of 
a full [hyper]sphere §d-i m any dimension d are characterized by the generic law (cos(#)) = 
exp(— (t/r)) where 6 denotes the generalized angular displacement of the walker^^. It was 
arguecP^that random walks on percolation clusters at threshold inscribed on [hyper] spheres 
would be characterized by relaxation of the form (cos(#)) = exp(— {t/r)" d ) with the same 
exponents (3d as in the corresponding Euclidean space. This was demonstrated numerically 
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for d = 3 to <*P^. A hypercube being topologically equivalent to a hypersphere, for random 
walks on a diluted hypercube at threshold one then expects stretched exponential relaxation 
with exponent (3 = 1/3. 

The diluted hypercube at threshold can alternatively be considered as a specific example 
of a sparse graph. Remarkably, analytic expressions for diffusion on general sparse graphs^^H 
derived from a quite different line of argument also lead to stretched exponential relaxation 
expressions with the same specific value (3 = 1/3 for the exponent. 

Here we present numerical data for random walks on the diluted hypercube at threshold 
up to dimension N = 28 which are consistent with these conclusions. We argue that the 
KWW relaxation observed phenomenologically in numerous complex systems just above 
their respective critical temperatures is not an artifact, but is the signature of a universal 
form of coarse grained configuration space morphology which precedes a glass transition. 



LAPLACE TRANSFORMS AND RANDOM NETWORKS 

Quite generally, any relaxation function q(t) can equivalently be characterized by its 
Laplacian, a relaxation mode density (or eigenvalue density) function p(s) defined by: 



POO 

q(t) = / p(s)e- st ds 
Jo 



with the normalization condition 



p(s)ds = 1 



(2) 



(3) 



In model systems it can be possible to establish analytically or numerically the distribution 
p(s) which can then be inverted to obtain q(t). The inverse Laplace transform of a numerical 
or experimental qx{t) to obtain p(s) is much more difficult unless q(t) is known to very 
high precision over a wide range of t. This is an ill-conditioned problem as different p(s) 
distributions can lead to almost indistinguishable q(t). 

PollarcP^ (see Berberan-Santos^) provided an exact inversion of the pure stretched ex- 
ponential relaxation function qx(t) [l] : 
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(4) 



For (3 < 1, pK,p(s) can be expressed in terms of elementary functions only for (3 = 1/2 251 ; in 
that case 



Pk,i/2(s) = [r^isrr^expi-l/AsT)] 



3/21 



(5) 



To a good approximation, for general (3 the large s (short time) limit takes the form px$ ~ 
s -{ l +P) anc i the small s (long time) limit the form p^^ ~ (s) expfs"^ 1 ^]. 

It should be kept in mind that at short times observed relaxation functions usually deviate 
from the "asymptotic" form. Also at very long times for finite sized systems the relaxation 
is controlled by the smallest non-zero value of s, s\. For time t > s^ 1 the relaxation will tend 
to a pure exponential, q(t) ~ exp[— tsi], but for large systems this condition corresponds to 
extremely long times and we will not consider it. What we are interested in is to establish 
the form of the relaxation in the regime where the mode distribution is no longer affected 
by short time effects and where p(s) can be considered continuous. 

RANDOM NETWORKS 

Random walks on the diluted iV-simplex or hypertetrahedron which is an Erdos-Renyi 
graph having dead ends and vertices with two connections, was studied theoretically by Bray 
and Rodgers^ using Replica theory. They showed that in this model the return function 
Pret{t), the probability that the walker will have returned to the origin after t steps, behaves 
like a stretched exponential with exponent 1/3. 

Samukhin et aP^ have made analytic studies of random walks and relaxation processes 
on uncorrelated Random Networks. They considered a stochastic process governed by the 
Laplacian operator occurring on a random graph with A* nodes, taking the limit as N* — > 
oo. They find that the determining parameter in this problem is the minimum degree q m of 
vertices (i.e. the minimum number of neighbors to any given vertex). For q m = 2, meaning 
that the network is "sparse" , the graph tends to a random Bethe lattice in which almost all 
finite subgraphs are trees, i.e., they contain almost no closed loops. In the present context 
the essential statement of Samukhin et aP^ is that when q m = 2 the mode density function 
Ps(s) for this very general model can be approximated by 

Ps(s) = s- 4/3 exp(-a/v^) (6) 

where 

-V? < 7 > 

with a similar expression for q m = 1 (graphs with dead ends). Then for a graph with N* 
vertices the asymptotics at t > In A* for the probability of return to the starting point at 
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time t during a random walk on the network (the "autocorrelator"^ will be 

Pr e t,5(t)~t ,? exp[-3(a/2) 2 / 3 t 1 / 3 ], (8) 

a stretched exponential having exponent 1/3, multiplied by a mildly time dependent prefac- 
tor (77 is small). This limit should be observable if the network size satisfies (lnA^*) 2//3 ^> 1. 

I. HYPERCUBE MODEL 

We have already addressed the hypercube problem numerically through Monte Carlo 
techniques^ and through the explicit solution of Master equation d 28 * 29 l In this paper we 
extend these results by investigating the time evolution for the autocorrelation function 
q(N,t), the return probability p ret (N,t), and the eigenvalue spectrum p(N,s) for diffusion 
on diluted hypercubes of dimension N near the critical occupation probability p c (N), for iV 
up to 28. 

Consider a hypercube (or n-cube) in [high] dimension N, Qn, with a fraction p of its 2 N 
vertices occupied at random. It is well establishe d 18 * 30 * 3 ^ that there is a critical threshold at 
p c {N) ~ 1/JV. For p > p c {N) the occupied vertices having one or more occupied vertices as 
neighbors make up a giant spanning cluster; for p < p c there exist only small clusters (each 
with less than N elements). By analogy with the equivalent situation in randomly occupied 
Euclidean space we will refer to p c as the "percolation" threshold. 

Gaunt and Brak® predict that the dependency of the critical site percolation concentra- 
tion p c on a hypercubic lattice of dimension d, Z d , or on a hypercube of dimension N, Qn, 
is given to order 4 by: 

Pc (a) = a+^a 2 + ^a 3 + 8 ^a 4 ... (9) 

where a(d) = 1/ (2d— 1) for the hypercubic lattice and cr(N) = 1/ (N—l) for the hypercube^. 
Although the terms in this expression are expected to be exact, the demonstration is not 
entirely rigorous^, and the series is obviously truncated. Grassberger^l tested the equation 
(JqJ) through large scale Monte Carlo simulations on Z d and verified that for d > 10 it 
represents the numerically determined p c {d) to within a small correction term. We will work 
with samples having vertex concentrations p(N) equal to the values p c (N) given by the 
truncated series equation ([9]). For different samples k the individual critical values p c (k) will 
in fact be distributed about the average valued 



For p > p c {N) we can define a random walk along edges on the giant cluster. Start at any 
vertex i on the giant cluster. Choose at random a vertex j on the hypercube, near neighbor 
to i. If the vertex j is also on the giant cluster and so accessible, move to j; otherwise the 
walker remains one time step longer at the vertex i. This evolution rule is chosen to mimic 
Monte Carlo simulations using Metropolis dynamics. 

We can compare the autocorrelation function q(N,t) obtained from this procedure, 
(q(N,t) is defined in Eq. ( Jl6| ) below), to the time dependent autocorrelation (Si(t).Si(0)) 
measured in thermodynamic models for systems of Ising spins SP and even to experimental 
magnetization decay results. From a theoretical point of view it is often more convenient 
to investigate the "return probability" p re t{t) that is basically the probability of finding the 



walker at the origin of the system after t steps (p ret (N,t) is defined in Eq. (18) below). For 
any network p re t(t) can be defined, while q(t) can be defined conveniently only on models 
such as the hypercube which have a suitable metric. 

The numerical data near criticality show that the long time relaxations of the autocorre- 
lation parameter q(N, t) and of the return probability p re t(N, t) are consistent with stretched 
exponentials having an exponent (3 = 1/3 over many orders of magnitude in time. 

ALGORITHM 

The time evolution of the entire probability distribution for the walker after t steps, n(t), 
can be described by a Master Equation. At t = the walker is localized on a single vertex 
i on the hypercube; the probability distribution then diffuses over the system at each time 
step following the equation: 



n i (t) = n i (t-i) + 



n»(t - i)w(j -M') - n,(t - i)w(i -> j) 



(10) 



where W (i — > j) represents the transition probability that is given by: 



if i vertex and j vertex are allowed 



W{i^j) = { N (11) 
otherwise 



The equation (10) can be rephrased as: 

n(t) = Fn(t-i) (12) 

6 



where F is the linear evolution operator. 

Since this process is Markovian we can diagonalize F; the smallest eigenvalue correspond- 
ing to the infinite time equilibrium limit (where all sites become equally populated) is 1. 
We can determine U and D satisfying: 

F = U T DU (13) 

where D is a diagonal matrix. For practical reasons it is convenient to diagonalize F so as 
to investigate the temporal evolution of the relevant quantities. We use: 

n(f) = F'n(o) = f/ T D*f/n(o) (14) 

We choose the initial condition as: 

n,(o) = 5 iio (15) 

where % is a vertex on the giant cluster. 

The value of the normalized autocorrelation function q(t) after time t for a given walk 
starting from i and arriving at i after time t can be defined by: 

* = feEE«^--) («>) 

\ io i I 

where dn(i,io) is the Hamming distance between vertex % and the initial state, No is the 
number of vertices on the giant cluster, q t=0 o for a given realization is given by: 

1 ^ N-2d H (i,io) h7) 
*=»-^Z, n (17) 

° iio 

and the averages are over different realizations of the diluted hypercube. 
We also calculated p ret defined by: 



We can show that: 



Pret(t) = — (}^ 1 -—) (19) 



This quantity is easier to calculate theoretically than q(t), but it is not useful to compare 
with results on model spin systems or experiments. We can write this equation in a more 
convenient form: 



where s = — In A and we excluded A = 1 eigenvalue. Another convenient form for investi- 
gating p ret is: 

poo 

Pret{t)= / dsp(s)e- ts (21) 
Jo 

where the density p is defined by: 




Our numerical workflow can be summarized as follows: 

1. generation of a diluted hypercube 

2. determination of the giant cluster 

3. determination of the eigenvalues and eigenvectors of F 

4. calculation of p(s), q{t) andp re t(t) 

The algorithm was implemented on Mathematica 8.0 and the simulations were performed 
on a Intel Xeon 2.27 Ghz with 24 Gbytes of Ram Memory. A single simulation for the N = 28 
cost 12 hours. The algorithm demands 24 Gbytes of memory for this case. 

Calculations were made with hypercubes of dimension N = 10, 12, 14, 16, 18, 20, 22, 24, 26 
and 28. All the calculations were performed at p c (N) values given by equation ([9j; this con- 
dition is important since it allows us to scale conveniently data for systems having different 
dimensions N. It is useful to be able to include data for smaller N in the global analysis as 
in these samples we deal with much smaller matrices which is simpler computationally. 

All vertices on the giant cluster were used as starting points, except for the largest 
systems N = 26 and 28 where we have not used all possible initial states i . For these sizes 
we approximated q(t) and p re t{t) by using only 1000 randomly chosen initial states for each 
realization. We have tested the accuracy of this approximation and we concluded that the 
error was very small (even for the smaller sizes). We studied 1000 different realizations of 
the hypercube for all sizes N except for iV = 28 when we have studied 100. 

II. NUMERICAL DATA 

On Figure ([I]) we represent a graphical representation of a diluted Q24 for this particular 
sample the graph is a tree showing the validity of the approximation proposed by^. 
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FIG. 1. A graphical representation of a diluted Q24 exactly at p c . The picture shows that the 
network presents no loops. 



The time evolution for the autocorrelation functions q(N,t) (16) is depicted in Figure 



[2] against log(t). On Figure [3] we show the equivalent results for the return probability 

Pret(N,t). 

In all cases we have fitted the long time part of the curves using the expression: 



f(t) = Aexp 



1/3 



(23) 



In Figures |4] and [5] we present the same results in a different manner so as to demonstrate 
the stretched exponential long time behavior. On the x axis the time scale is normalized with 
x(t) = (t/r(A^)) 1 / 3 and on the y axis the measured q(N,t) or p ret (N,t) are normalized so 
y(N,t) = ]n(q(N,t)/A q (N)) and y(N,t) = ln(p ret (N,t)/A ret (N)) respectively. In these plots 
a stretched exponential with exponent 1/3 is a straight line as observed; we have chosen 
the normalization factors r(N) and A q (N), A ret (N) so that data for different hypercube 
dimensions N collapse. This form of plot allows one to distinguish clearly between the short 
time regime and the stretched exponential regime; the latter can be seen to extend over 
a wide time range until measurements are limited by the statistical noise. The effective 
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FIG. 3. The decay of the return probability logp re t(N,t), Eqn.(18), against logt for N from 10 
to 28. 
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FIG. 4. The decay of the normalized autocorrelation function ln(q(N, t)/A q (N)) against (i/r) 1 / 3 . 
For stretched exponentials with exponent (3 = 1/3 in the long time regime the data should lie on 
a straight line in this form of plot, as observed. 

exponent (3 = 1/3 is independent of TV to within the statistics. 

On Figure [6] we show the size dependence on the t(N) time scale parameter from the fits 
of the autocorrelation q(t,N) and the return probability p ret (t,N) data. The data can be 
fitted by fitted by 



with the fit parameters 7 = 0.24 ± 0.1 and B = 1.5 ± 0.1 for autocorrelation function, 
7 = 0.24 ± 0.05 and B = 1.7 ± 0.2 for the return probability. The values of the time scaling 
parameters t(N) for the two different observables are identical within the precision of the 
measurements. 

The most fundamental way to understand the system dynamics is through investigating 
the eigenvalue spectra; the stretched exponential long time behavior depends exclusively 
on the density of the eigenvalues above the smallest eigenvalue, in the region where the 
distribution for a finite size sample can still be considered to be continuous. A given spectrum 
leads unambiguously to a unique relaxation function, while it is much more difficult to 



t(n) = mo 77V 



(24) 
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FIG. 5. The decay of the normalized return probability log(p re t(N,t)/A re t(N)) against 
(t/r(A r )) 1 / 3 . For stretched exponentials with exponent (3 = 1/3 in the long time regime the 
data should lie on a straight line on this form of plot, as observed. 

determine the precise form of a mode spectrum from a relaxation function. 

On Figure [7] we compare the mode density p(s) obtained through the present simulations 
with the theoretical expressions. All the numerical results were obtained using 1000 different 
realizations of the diluted hypercube at each dimension N. Unfortunately in practice the 
calculations of p(s) are numerically demanding because of strong sample to sample fluctua- 
tions. The spectra were first binned in the form of histograms. We defined a cut-off \ m in{N) 
or equivalently s max (N) = — \n\ min (N) to eliminate the short time effects and selected the 
eigenvalues on the interval s G (0,s max (N)). We choose s max = 2/t(N) for all dimensions. 
We divided this interval in bins equally spaced on a logarithmic scale and then calculated 
the densities for each interval, normalizing the frequencies by the length of the intervals. 

The continuous curves were calculated from the expression Q for pK,i/z{^) arid from the 



approximate analytic expression (J6J) for ps(A) using t(N) estimated from equation (24). To 
compare with simulation results we normalized p(A) functions using: 
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FIG. 6. The dependence of the time scale t(N) with dimension for the return probability p re t(N, t) 
(in red) and autocorrelation q(N,t) (in blue). 

C- 1 = / p(s)ds (25) 
Jo 

and 

p'(s) = Cp(s) (26) 

Over the ranges for which reliable data points have been obtained the measured mode 
spectrum densities p(N, s) closely resemble the corresponding parts of the calculated spectra 
from the Laplace transform Pk,i/^{ s )i @ or the analytic ps(s) spectrum (which are in 
fact very similar to each other). The numerical spectra for the hypercube model are indeed 
consistent with the mode density spectral form derived analytically for the more general 
random network modeP^. 

DISCUSSION AND CONCLUSIONS 

We have studied numerically relaxation through random walks along near neighbor edges 
on the giant cluster of vertices in randomly diluted hypercubes of dimensions up to iV = 28 
near the percolation threshold for the cluster. The data show clearly that at the percolation 
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FIG. 7. Spectral density data p(N,s) from the hypercube evaluations together with the exact 
Laplace transforms Pi/s(s) Q for stretched exponentials with /3 = 1/3 and t(N) values equal to 
the numerical estimates 24 (dashed lines), and the analytic sparse network expression (j6|P (full 
lines). The values are normalized (see text). 



threshold concentration p c (N), the relaxation mode spectrum, the time dependence of the 
autocorrelation q(N,t), and the return probability p ret (N,t), are all consistent with asymp- 
totic stretched exponential relaxation exp[— (t/r^N)) 13 } having exponent (3 = 1/3. The time 



scale t(N) increases exponentially with dimension N, Eqn. (24). The observed eigenvalue 
spectra demonstrate that the dynamical q(N,t) behavior previously obtained from Monte 
Carlo simulations and from numerical solutions of the master equatioEpZES ^ oes no t repre- 
sent a crossover between different exponential regimes, but that it is the consequence of a 
specific wide eigenvalue spectrum. 

A final long time crossover to a pure exponential (which would correspond to a regime 
where the effective relaxation mode spectrum is reduced to a gap between the ground state 
and the lowest mode) is not visible in the data. 

This diluted hypercube model at threshold can be considered as the limiting high di- 
mensional case of percolation on sphere-like spaces. Alternatively it can be considered as 
a specific explicit example of a generic sparse random network. The observed stretched 
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exponential behavior with exponent (3 = 1/3 on the dilute hypercube at the percolation 
threshold is consistent with the predictions of the sphere-like percolation approach^ and 
with studies of random walks on sparse random network d 23 * 24 ! , where the same stretched 
exponential relaxation with the same exponent — 1/3 has been derived analytically. 

For a physical system, configuration space can be imagined as a very high dimensional 
graph. The system's dynamics is equivalent to a random walk of the point representing the 
instantaneous state of the system among those vertices of the graph which are thermody- 
namically accessible. We suggest that when the stretched exponential exp[— (t/r) 1 / 3 ] form 
of limiting relaxation with diverging r is observed numerically or experimentally for the au- 
tocorrelation function relaxation q(t) in complex physical systems (which is often the case, 
see for instance ^ 35 * 36 l)it is the signature of a configuration space tending to a percolation 
threshold and having a sparse random network topology. 
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